duration = c(seq.Date(peace_talks$start[1], peace_talks$end[1], by = "day"),
seq.Date(peace_talks$start[2], peace_talks$end[2], by = "day"),
seq.Date(peace_talks$start[3], peace_talks$end[3], by = "day"),
seq.Date(peace_talks$start[4], peace_talks$end[4], by = "day"),
seq.Date(peace_talks$start[5], peace_talks$end[5], by = "day"),
seq.Date(peace_talks$start[6], peace_talks$end[6], by = "day"))
length(duration)
tishreen1 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[1] & date < peace_talks$end[1], 1,
ifelse(date %in% duration, NA, 0)))
tishreen2 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[2] & date < peace_talks$end[2], 1,
ifelse(date %in% duration, NA, 0)))
tishreen3 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[3] & date < peace_talks$end[3], 1,
ifelse(date %in% duration, NA, 0)))
tishreen4 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[4] & date < peace_talks$end[4], 1,
ifelse(date %in% duration, NA, 0)))
tishreen5 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[5] & date < peace_talks$end[5], 1,
ifelse(date %in% duration, NA, 0)))
tishreen6 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[6] & date < peace_talks$end[6], 1,
ifelse(date %in% duration, NA, 0)))
library(estimatr)
library(texreg)
mod1 = lm_robust(isr_pct ~ talks + year, data = tishreen1)
mod2 = lm_robust(isr_pct ~ talks + year, data = tishreen2)
mod3 = lm_robust(isr_pct ~ talks + year, data = tishreen3)
mod4 = lm_robust(isr_pct ~ talks + year, data = tishreen4)
mod5 = lm_robust(isr_pct ~ talks + year, data = tishreen5)
mod6 = lm_robust(isr_pct ~ talks + year, data = tishreen6)
texreg(list(mod1, mod2, mod3, mod4, mod5, mod6), custom.model.names = c("1989/12-1991/10", "1992/07-1993/09", "1994/04-1994/12",
"1995/05-1995/11", "1995/12-1996/02", "1999/12-2000/01"),
omit.coef = "year", custom.coef.names = c("Intercept", "Peace Talks"), float.pos = "H",
caption = "Coverage of Israel in Tishreen during peace talks with year fixed effects", label = "tab:peace_talks_appendix",
fontsize = "footnotesize", include.ci = F, include.rsquared = F, include.adjrs = F, include.nobs = F,
reorder.coef = c(2, 1), include.rmse = F, file = "tabs/appendix/peace_talks_linear.tex")
# permutations ------------------------------------------------------------
mod1 = lm_robust(isr_pct ~ talks + year, data = tishreen1)
mod2 = lm_robust(isr_pct ~ talks + year, data = tishreen2)
mod3 = lm_robust(isr_pct ~ talks + year, data = tishreen3)
mod4 = lm_robust(isr_pct ~ talks + year, data = tishreen4)
mod5 = lm_robust(isr_pct ~ talks + year, data = tishreen5)
mod6 = lm_robust(isr_pct ~ talks + year, data = tishreen6)
beta1 = coef(lm(isr_pct ~ talks, data = tishreen1))[2]
beta2 = coef(lm(isr_pct ~ talks, data = tishreen2))[2]
beta3 = coef(lm(isr_pct ~ talks, data = tishreen3))[2]
beta4 = coef(lm(isr_pct ~ talks, data = tishreen4))[2]
beta5 = coef(lm(isr_pct ~ talks, data = tishreen5))[2]
beta6 = coef(lm(isr_pct ~ talks, data = tishreen6))[2]
nsims = 1e4
store_betas = data.frame(coef1 = numeric(nsims), coef2 = numeric(nsims), coef3 = numeric(nsims),
coef4 = numeric(nsims), coef5 = numeric(nsims), coef6 = numeric(nsims))
minmonth = peace_talks$start[1]
maxmonth = peace_talks$start[6]
dates = seq.Date(minmonth, maxmonth, by = "day")
s1 = length(seq.Date(peace_talks$start[1], peace_talks$end[1], by = "day")) - 1
s2 = length(seq.Date(peace_talks$start[2], peace_talks$end[2], by = "day")) - 1
s3 = length(seq.Date(peace_talks$start[3], peace_talks$end[3], by = "day")) - 1
s4 = length(seq.Date(peace_talks$start[4], peace_talks$end[4], by = "day")) - 1
s5 = length(seq.Date(peace_talks$start[5], peace_talks$end[5], by = "day")) - 1
s6 = length(seq.Date(peace_talks$start[6], peace_talks$end[6], by = "day")) - 1
nsim = 1e2
set.seed(32453245)
# do the simulations
for (i in 1:nsims){
# select treatment dates to kick in
start = sample(dates, 1)
sampled1 = seq.Date(start, start + s1, by = "day")
sampled2 = seq.Date(start, start + s2, by = "day")
sampled3 = seq.Date(start, start + s3, by = "day")
sampled4 = seq.Date(start, start + s4, by = "day")
sampled5 = seq.Date(start, start + s5, by = "day")
sampled6 = seq.Date(start, start + s6, by = "day")
d1 = tishreen %>% mutate(treated = ifelse(date %in% sampled1, 1, ifelse(date %in% duration, NA, 0)))
d2 = tishreen %>% mutate(treated = ifelse(date %in% sampled2, 1, ifelse(date %in% duration, NA, 0)))
d3 = tishreen %>% mutate(treated = ifelse(date %in% sampled3, 1, ifelse(date %in% duration, NA, 0)))
d4 = tishreen %>% mutate(treated = ifelse(date %in% sampled4, 1, ifelse(date %in% duration, NA, 0)))
d5 = tishreen %>% mutate(treated = ifelse(date %in% sampled5, 1, ifelse(date %in% duration, NA, 0)))
d6 = tishreen %>% mutate(treated = ifelse(date %in% sampled6, 1, ifelse(date %in% duration, NA, 0)))
# run regression
mod1 = lm(isr_pct ~ treated, data = d1)
mod2 = lm(isr_pct ~ treated, data = d2)
mod3 = lm(isr_pct ~ treated, data = d3)
mod4 = lm(isr_pct ~ treated, data = d4)
mod5 = lm(isr_pct ~ treated, data = d5)
mod6 = lm(isr_pct ~ treated, data = d6)
store_betas$coef1[i] = coef(mod1)["treated"]
store_betas$coef2[i] = coef(mod2)["treated"]
store_betas$coef3[i] = coef(mod3)["treated"]
store_betas$coef4[i] = coef(mod4)["treated"]
store_betas$coef5[i] = coef(mod5)["treated"]
store_betas$coef6[i] = coef(mod6)["treated"]
}
# get one-sided p-value
pval1 = mean(store_betas$coef1 < beta1, na.rm = T)
pval2 = mean(store_betas$coef2 < beta2, na.rm = T)
pval3 = mean(store_betas$coef3 < beta3, na.rm = T)
pval4 = mean(store_betas$coef4 < beta4, na.rm = T)
pval5 = mean(store_betas$coef5 < beta5, na.rm = T)
pval6 = mean(store_betas$coef6 < beta6, na.rm = T)
pval1
pval2
pval3
pval4
pval5
pval6
# tishreen plot results ------------------------------------------------------------
p1 = ggplot(store_betas) +
aes(x = coef1) +
geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) +
geom_segment(x = beta1, xend = beta1, y = 0, yend = .07, size = 1) +
annotate(geom = "text", x = beta1, y = .08, label = paste0("1989-1991\nOne-sided\np = ", round(pval1, 3))) +
labs(x = "Placebo coefficient estimate", y = "Density") +
scale_x_continuous(breaks = seq(-.1, .1, .05)) +
theme_minimal() +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=.85),
panel.spacing = unit(0, "in"))
p1
#This code produces Table 15 in appendix F and figures 20, 21, 22, 23, 24, and 25 in Appendix G
if(!dir.exists("tabs")){dir.create("tabs")}
if(!dir.exists("tabs/appendix")){dir.create("tabs/appendix")}
if(!dir.exists("figs/appendix")){dir.create("figs/appendix")}
rm(list = ls())
library(lubridate)
library(foreign)
library(ggthemes)
library(stats)
library(lmtest)
library(sandwich)
library(quanteda)
library(stm)
library(stats)
library(ggpubr)
library(tidyverse)
library(tidytext)
gc()
# Load tishreen -----------------------------------------------------------
tishreen = read.csv("tishreen.csv")
tishreen = tishreen %>%
as_tibble() %>%
mutate(date = ymd(date),
week = floor_date(date, unit = "week"),
isr_pct = ofhlisrael/totalofhl) %>%
group_by(week) %>%
mutate(isr_avg = sum(ofhlisrael)/sum(totalofhl))
#Generate missing data
df = tibble(date = seq(ymd("1987-01-03"), ymd("2002-12-31"), by = "1 day"))
tishreen = left_join(df, tishreen, "date")
tishreen$year = as.factor(floor_date(tishreen$date, unit = "year"))
peace_talks = tibble(start = c(ymd("1989-12-01"), ymd("1992-07-01"), ymd("1994-04-30"), ymd("1995-05-24"), ymd("1995-12-11"), ymd("1999-12-08")),
end = c(ymd("1991-10-30"), ymd("1993-09-13"), ymd("1994-12-23"), ymd("1995-11-04"), ymd("1996-02-25"), ymd("2000-01-11")))
duration = c(seq.Date(peace_talks$start[1], peace_talks$end[1], by = "day"),
seq.Date(peace_talks$start[2], peace_talks$end[2], by = "day"),
seq.Date(peace_talks$start[3], peace_talks$end[3], by = "day"),
seq.Date(peace_talks$start[4], peace_talks$end[4], by = "day"),
seq.Date(peace_talks$start[5], peace_talks$end[5], by = "day"),
seq.Date(peace_talks$start[6], peace_talks$end[6], by = "day"))
length(duration)
tishreen1 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[1] & date < peace_talks$end[1], 1,
ifelse(date %in% duration, NA, 0)))
tishreen2 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[2] & date < peace_talks$end[2], 1,
ifelse(date %in% duration, NA, 0)))
tishreen3 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[3] & date < peace_talks$end[3], 1,
ifelse(date %in% duration, NA, 0)))
tishreen4 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[4] & date < peace_talks$end[4], 1,
ifelse(date %in% duration, NA, 0)))
tishreen5 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[5] & date < peace_talks$end[5], 1,
ifelse(date %in% duration, NA, 0)))
tishreen6 = tishreen %>% mutate(talks = ifelse(date > peace_talks$start[6] & date < peace_talks$end[6], 1,
ifelse(date %in% duration, NA, 0)))
library(estimatr)
library(texreg)
mod1 = lm_robust(isr_pct ~ talks + year, data = tishreen1)
mod2 = lm_robust(isr_pct ~ talks + year, data = tishreen2)
mod3 = lm_robust(isr_pct ~ talks + year, data = tishreen3)
mod4 = lm_robust(isr_pct ~ talks + year, data = tishreen4)
mod5 = lm_robust(isr_pct ~ talks + year, data = tishreen5)
mod6 = lm_robust(isr_pct ~ talks + year, data = tishreen6)
texreg(list(mod1, mod2, mod3, mod4, mod5, mod6), custom.model.names = c("1989/12-1991/10", "1992/07-1993/09", "1994/04-1994/12",
"1995/05-1995/11", "1995/12-1996/02", "1999/12-2000/01"),
omit.coef = "year", custom.coef.names = c("Intercept", "Peace Talks"), float.pos = "H",
caption = "Coverage of Israel in Tishreen during peace talks with year fixed effects", label = "tab:peace_talks_appendix",
fontsize = "footnotesize", include.ci = F, include.rsquared = F, include.adjrs = F, include.nobs = F,
reorder.coef = c(2, 1), include.rmse = F, file = "tabs/appendix/peace_talks_linear.tex")
# permutations ------------------------------------------------------------
mod1 = lm_robust(isr_pct ~ talks + year, data = tishreen1)
mod2 = lm_robust(isr_pct ~ talks + year, data = tishreen2)
mod3 = lm_robust(isr_pct ~ talks + year, data = tishreen3)
mod4 = lm_robust(isr_pct ~ talks + year, data = tishreen4)
mod5 = lm_robust(isr_pct ~ talks + year, data = tishreen5)
mod6 = lm_robust(isr_pct ~ talks + year, data = tishreen6)
beta1 = coef(lm(isr_pct ~ talks, data = tishreen1))[2]
beta2 = coef(lm(isr_pct ~ talks, data = tishreen2))[2]
beta3 = coef(lm(isr_pct ~ talks, data = tishreen3))[2]
beta4 = coef(lm(isr_pct ~ talks, data = tishreen4))[2]
beta5 = coef(lm(isr_pct ~ talks, data = tishreen5))[2]
beta6 = coef(lm(isr_pct ~ talks, data = tishreen6))[2]
nsims = 1e2
store_betas = data.frame(coef1 = numeric(nsims), coef2 = numeric(nsims), coef3 = numeric(nsims),
coef4 = numeric(nsims), coef5 = numeric(nsims), coef6 = numeric(nsims))
minmonth = peace_talks$start[1]
maxmonth = peace_talks$start[6]
dates = seq.Date(minmonth, maxmonth, by = "day")
s1 = length(seq.Date(peace_talks$start[1], peace_talks$end[1], by = "day")) - 1
s2 = length(seq.Date(peace_talks$start[2], peace_talks$end[2], by = "day")) - 1
s3 = length(seq.Date(peace_talks$start[3], peace_talks$end[3], by = "day")) - 1
s4 = length(seq.Date(peace_talks$start[4], peace_talks$end[4], by = "day")) - 1
s5 = length(seq.Date(peace_talks$start[5], peace_talks$end[5], by = "day")) - 1
s6 = length(seq.Date(peace_talks$start[6], peace_talks$end[6], by = "day")) - 1
set.seed(32453245)
# do the simulations
for (i in 1:nsims){
# select treatment dates to kick in
start = sample(dates, 1)
sampled1 = seq.Date(start, start + s1, by = "day")
sampled2 = seq.Date(start, start + s2, by = "day")
sampled3 = seq.Date(start, start + s3, by = "day")
sampled4 = seq.Date(start, start + s4, by = "day")
sampled5 = seq.Date(start, start + s5, by = "day")
sampled6 = seq.Date(start, start + s6, by = "day")
d1 = tishreen %>% mutate(treated = ifelse(date %in% sampled1, 1, ifelse(date %in% duration, NA, 0)))
d2 = tishreen %>% mutate(treated = ifelse(date %in% sampled2, 1, ifelse(date %in% duration, NA, 0)))
d3 = tishreen %>% mutate(treated = ifelse(date %in% sampled3, 1, ifelse(date %in% duration, NA, 0)))
d4 = tishreen %>% mutate(treated = ifelse(date %in% sampled4, 1, ifelse(date %in% duration, NA, 0)))
d5 = tishreen %>% mutate(treated = ifelse(date %in% sampled5, 1, ifelse(date %in% duration, NA, 0)))
d6 = tishreen %>% mutate(treated = ifelse(date %in% sampled6, 1, ifelse(date %in% duration, NA, 0)))
# run regression
mod1 = lm(isr_pct ~ treated, data = d1)
mod2 = lm(isr_pct ~ treated, data = d2)
mod3 = lm(isr_pct ~ treated, data = d3)
mod4 = lm(isr_pct ~ treated, data = d4)
mod5 = lm(isr_pct ~ treated, data = d5)
mod6 = lm(isr_pct ~ treated, data = d6)
store_betas$coef1[i] = coef(mod1)["treated"]
store_betas$coef2[i] = coef(mod2)["treated"]
store_betas$coef3[i] = coef(mod3)["treated"]
store_betas$coef4[i] = coef(mod4)["treated"]
store_betas$coef5[i] = coef(mod5)["treated"]
store_betas$coef6[i] = coef(mod6)["treated"]
}
# get one-sided p-value
pval1 = mean(store_betas$coef1 < beta1, na.rm = T)
pval2 = mean(store_betas$coef2 < beta2, na.rm = T)
pval3 = mean(store_betas$coef3 < beta3, na.rm = T)
pval4 = mean(store_betas$coef4 < beta4, na.rm = T)
pval5 = mean(store_betas$coef5 < beta5, na.rm = T)
pval6 = mean(store_betas$coef6 < beta6, na.rm = T)
pval1
pval2
pval3
pval4
pval5
pval6
# tishreen plot results ------------------------------------------------------------
p1 = ggplot(store_betas) +
aes(x = coef1) +
geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) +
geom_segment(x = beta1, xend = beta1, y = 0, yend = .07, size = 1) +
annotate(geom = "text", x = beta1, y = .08, label = paste0("1989-1991\nOne-sided\np = ", round(pval1, 3))) +
labs(x = "Placebo coefficient estimate", y = "Density") +
scale_x_continuous(breaks = seq(-.1, .1, .05)) +
theme_minimal() +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=.85),
panel.spacing = unit(0, "in"))
p1
ggsave("figs/appendix/tishreen_placebo1.pdf", width=11, height = 7)
p2 = ggplot(store_betas) +
aes(x = coef2) +
geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) +
geom_segment(x = beta2, xend = beta2, y = 0, yend = .07, size = 1) +
annotate(geom = "text", x = beta2 - 0.02, y = .074, label = paste0("1992-1993\nOne-sided\np = ", round(pval2, 3))) +
labs(x = "Placebo coefficient estimate", y = "Density") +
scale_x_continuous(breaks = seq(-.1, .1, .05)) +
theme_minimal() +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=.85),
panel.spacing = unit(0, "in"))
ggsave("figs/appendix/tishreen_placebo2.pdf", width=11, height = 7)
p3 = ggplot(store_betas) +
aes(x = coef3) +
geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) +
geom_segment(x = beta3, xend = beta3, y = 0, yend = .07, size = 1) +
annotate(geom = "text", x = beta3 - 0.02, y = .074, label = paste0("1994-1994\nOne-sided\np = ", round(pval3, 3))) +
labs(x = "Placebo coefficient estimate", y = "Density") +
scale_x_continuous(breaks = seq(-.1, .1, .05)) +
theme_minimal() +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=.85),
panel.spacing = unit(0, "in"))
p3
ggsave("figs/appendix/tishreen_placebo3.pdf", width=11, height = 7)
p4 = ggplot(store_betas) +
aes(x = coef4) +
geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) +
geom_segment(x = beta4, xend = beta4, y = 0, yend = .07, size = 1) +
annotate(geom = "text", x = beta4 - 0.02, y = .074, label = paste0("1995-1995\nOne-sided\np = ", round(pval4, 3))) +
labs(x = "Placebo coefficient estimate", y = "Density") +
scale_x_continuous(breaks = seq(-.1, .1, .05)) +
theme_minimal() +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=.85),
panel.spacing = unit(0, "in"))
p4
ggsave("figs/appendix/tishreen_placebo4.pdf", width=11, height = 7)
p5 = ggplot(store_betas) +
aes(x = coef5) +
geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) +
geom_segment(x = beta5, xend = beta5, y = 0, yend = .07, size = 1) +
annotate(geom = "text", x = beta5 - 0.03, y = .074, label = paste0("1995-1996\nOne-sided\np = ", round(pval5, 3))) +
labs(x = "Placebo coefficient estimate", y = "Density") +
scale_x_continuous(breaks = seq(-.1, .1, .05)) +
theme_minimal() +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=.85),
panel.spacing = unit(0, "in"))
p5
ggsave("figs/appendix/tishreen_placebo5.pdf", width=11, height = 7)
p6 = ggplot(store_betas) +
aes(x = coef6) +
geom_histogram(fill = "#A0A0A0", aes(y = ..ncount../sum(..ncount..)), bins = 40) +
geom_segment(x = beta6, xend = beta6, y = 0, yend = .07, size = 1) +
annotate(geom = "text", x = beta6 - 0.03, y = .074, label = paste0("1999-2000\nOne-sided\np = ", round(pval6, 3))) +
labs(x = "Placebo coefficient estimate", y = "Density") +
scale_x_continuous(breaks = seq(-.1, .1, .05)) +
theme_minimal() +
theme(panel.grid = element_blank(),
panel.border = element_rect(colour = "black", fill=NA, size=.85),
panel.spacing = unit(0, "in"))
p6
ggsave("figs/appendix/tishreen_placebo6.pdf", width=11, height = 7)
#This code produces Table 16 in Appendix H
if(!dir.exists("tabs")){dir.create("tabs")}
if(!dir.exists("tabs/appendix")){dir.create("tabs/appendix")}
rm(list = ls())
library(lubridate)
library(foreign)
library(ggthemes)
library(stats)
library(lmtest)
library(sandwich)
library(quanteda)
library(stm)
library(stats)
library(ggpubr)
library(tidyverse)
library(tidytext)
# tishreen ----------------------------------------------------------------
tishreen = read.csv("tishreen.csv")
tishreen = tishreen %>%
as_tibble() %>%
mutate(date = ymd(date),
week = floor_date(date, unit = "week"),
isr_pct = ofhlisrael/totalofhl) %>%
group_by(week) %>%
mutate(isr_avg = sum(ofhlisrael)/sum(totalofhl))
#Generate missing data
df = tibble(date = seq(ymd("1987-01-03"), ymd("2002-12-31"), by = "1 day"))
tishreen = left_join(df, tishreen, "date")
# weekdays tishreen ----------------------------------------------------------------
library(estimatr)
library(texreg)
mod1 = lm_robust(ofhlisrael ~ ofpicturesofbasharorhafez, data = tishreen)
mod2 = lm_robust(ofhlisrael ~ ofpicturesofbasharorhafez + as.factor(floor_date(date, "year")), data = tishreen)
texreg(list(mod1, mod2), omit.coef = "factor", custom.coef.names = c("Intercept", "Assad Pictures"),
caption = "Relationship between coverage of Israel and pictures of Assad on \\textit{Tishreen} front page",
label = "tab:assad_pics", include.ci = F, include.rsquared = F, include.adjrs = F,
include.nobs = F, reorder.coef = c(2, 1), include.rmse = F,
float.pos = "H", custom.gof.rows = list("Year FE" = c("No", "Yes")))
#This code produces Table 16 in Appendix H
if(!dir.exists("tabs")){dir.create("tabs")}
if(!dir.exists("tabs/appendix")){dir.create("tabs/appendix")}
rm(list = ls())
library(lubridate)
library(foreign)
library(ggthemes)
library(stats)
library(lmtest)
library(sandwich)
library(quanteda)
library(stm)
library(stats)
library(ggpubr)
library(tidyverse)
library(tidytext)
# tishreen ----------------------------------------------------------------
tishreen = read.csv("tishreen.csv")
tishreen = tishreen %>%
as_tibble() %>%
mutate(date = ymd(date),
week = floor_date(date, unit = "week"),
isr_pct = ofhlisrael/totalofhl) %>%
group_by(week) %>%
mutate(isr_avg = sum(ofhlisrael)/sum(totalofhl))
#Generate missing data
df = tibble(date = seq(ymd("1987-01-03"), ymd("2002-12-31"), by = "1 day"))
tishreen = left_join(df, tishreen, "date")
# weekdays tishreen ----------------------------------------------------------------
library(estimatr)
library(texreg)
mod1 = lm_robust(ofhlisrael ~ ofpicturesofbasharorhafez, data = tishreen)
mod2 = lm_robust(ofhlisrael ~ ofpicturesofbasharorhafez + as.factor(floor_date(date, "year")), data = tishreen)
texreg(list(mod1, mod2), omit.coef = "factor", custom.coef.names = c("Intercept", "Assad Pictures"),
caption = "Relationship between coverage of Israel and pictures of Assad on \\textit{Tishreen} front page",
label = "tab:assad_pics", include.ci = F, include.rsquared = F, include.adjrs = F,
include.nobs = F, reorder.coef = c(2, 1), include.rmse = F,
float.pos = "H", custom.gof.rows = list("Year FE" = c("No", "Yes")), file = "tabs/appendix/assad_pics")
files = list.files(pattern="[.]R$")
files
files = files[!grepl("*source_all.R", files)]
files
for(f in files){
message("starting file ", f)
sink(paste0(f, ".txt"))
source(f, echo=T)
sink()
}
if(!dir.exists("figs")){dir.create("figs")}
if(!dir.exists("tabs")){dir.create("tabs")}
library(tidyverse)
library(arabicStemR)
library(tidytext)
library(topicmodels)
library(lubridate)
library(ggplot2)
library(quanteda)
library(stm)
library(stats)
library(ggthemes)
library(ggpubr)
rm(list = ls())
# Reading data ------------------------------------------------------------
load("k40.RData")
beta = tidy(topic_model)
gamma = tidy(topic_model, matrix = "gamma")
meta = meta %>% mutate(document = article) %>% select(- article)
gamma = left_join(gamma, meta, by = "document")
gc()
# Naming topics -----------------------------------------------------------
topic_avg = gamma %>%
group_by(topic) %>%
summarise(avg = mean(gamma)) %>%
arrange(-avg)
terms = labelTopics(topic_model, n = 30)
topic_terms = terms$prob %>% as_tibble %>% mutate(topic = row_number()) %>% unite("terms", V1:V30, sep = ",") %>% select(topic, terms)
df = topic_terms %>% left_join(topic_avg) %>%
mutate(topic_name = c("communication", "assad", "isis", "temperature", "elections",
"weather", "foreign fighters", "united states", "announcements", "conspiracies and plots",
"election winners", "legislation", "terrorist attacks", "speeches", "victims",
"gulf", "accident", "media", "europe", "yemen",
"culture and sports", "capitulation", "economy", "diplomacy", "religion",
"transport", "natural resources", "national unity", "lebanon", "finance",
"bureaucracy", "education", "fighting terrorism", "Israel and Palestine", "ba'th party",
"Iraq", "Turkey", "russia", "Iran", "regional politics")) %>%
mutate(topic_name = str_to_title(topic_name)) %>%
select(topic_name, terms, avg) %>%
rename("Topic Label" = "topic_name", "High Proabbility Terms" = "terms", "Expected Proportion" = "avg")
df = df %>% mutate(`High Proabbility Terms` = c("add, call, say, information, electronic, details",
"president, mister, Assad, Bashar, Arab, Syria",
"terrorism, organization, Syria, ISIS, state, armed",
"[regions in Syria], degree, temperature, increase, averages, high, low",
"election, constitution, people, presidency, council, candidate",
"[regions in Syria], agriculture, reach, rain, season",
"Newspaper, British, Tunisian, French, intelligence",
"America, United, States, president, Washington, military",
"first, public, SANA, yesterday, thawrah, Damascus",
"people, resistance, Zionist, enemy, conpiracy, interference",
"[winners' names], regional, doctor, MP, council, committee",
"article, law, decree, ruling, number, court",
"armed, group, terrorist, elements, fire",
"said, politics, states, should, change, discuss, people",
"aid, humanitarian, crisis, support, hungry, camp",
"Saudi, Egypt, Qatar, regime, gulf, opposition, kingdom",
"bombing, crime, children, civilian, assault",
"media, ask, truth, channel, report, events, news",
"Europe, union, France, nations, Germany, nuclear",
"Yemen, enemy, Saudi, army, raid, bombing",
"culture, Syria, world, exhibition, modern, sport",
"Aleppo, city, people, army, terrorism, reconciliation, Homs",
"economy, Syria, sector, trade, invest, tourism",
"relation, visit, state, delegation, meeting, foreign",
"religion, Islam, Patriarch, peace, Christian, Muslim",
"Jordan, sea, fly, border, transport, airport",
"oil, water, project, electricity, station, gas",
"nation, Syria, martyrs, army, sons, people, land, victory",
"Lebanon, resistance, Israel, president, nation, enemy",
"money, public, firm, Lira, bank, dollar",
"council, project, work, public, manage, minister",
"university, education, student, study, higher, public",
"terror, source, army, armed, organization, destruction, kill",
"Israel, Palestine, occupation, strip (Gaza), Jerusalem, settler",
"nation (watan; qawm), party, Arab, leadership, people, Ba'th",
"Iraq, Baghdad, kill, America, injure, city",
"Turkey, Erdogan, party, government, regime, development",
"Russia(n), Moscow, united, Foreign, Lavrov",
"Iran, state, region, Islam, Tehran, nuclear",
"state, Arab, Israel, peace, decision, conference"))
df %>% arrange(desc(`Expected Proportion`)) %>%
xtable::xtable(caption = "Topic labels, high probability terms for each topic, and their expected proportion",
label = "tab:topics") %>%
xtable::print.xtable(include.rownames = F, file = "tabs/topics.tex")
# Topics proportions ------------------------------------------------------
library(ggplot2)
library(ggthemes)
df %>%
arrange(desc(`Expected Proportion`)) %>%
mutate(`Topic Label` = factor(`Topic Label`, levels = `Topic Label`),
`Topic Label` = factor(`Topic Label`, levels = rev(levels(`Topic Label`)))) %>%
ggplot(aes(x = `Topic Label`, y = `Expected Proportion`)) +
geom_col(width = .1) +
coord_flip() +
theme_few() +
labs(y = "Expected Topic Proportions", x = "")
ggsave("figs/topic_proportions.pdf", height = 6, width = 6)
